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Across many fields, a problem of interest is to predict the transition rates between nodes of a 
network, given limited stationary state and dynamical information. We give a solution using the 
principle of Maximum Caliber. We find the transition rate matrix by maximizing the path entropy 
of a random walker on the network constrained to reproducing a stationary distribution and a few 
dynamical averages. A main finding here is that when constrained only by the mean jump rate, 
the rate matrix is given by a square-root dependence of the rate, Wab oc ^/pb/pa, on pa and pb, the 
stationary state populations at nodes a and b. We give two examples of our approach. First, we 
show that this method correctly predicts the correlated rates in a biochemical network of two genes, 
where we know the exact results from prior simulation. Second, we show that it correctly predicts 
rates of peptide conformational transitions, when compared to molecular dynamics simulations. 
This method can be used to infer large numbers of rates on known networks where smaller numbers 
of steady-state node populations are known. 

INTRODUCTION 

We are interested in an inference problem in network 
science. Given the topology of a network and station¬ 
ary populations at the nodes, what is the best model 
that can infer the rates of the dynamical flows along the 
edges? Here are examples. First, consider a spin model 
with a known stationary distribution, for example, those 
used in neuroscience o, protein evolution m, or col¬ 
loidal sciences ©■ It is of great interest to infer the best 
dynamical process that is consistent with a given rate 
of spin flip. Second, in systems biology, we often know 
the topology of a network of metabolites, or proteins or 
regulatory elements. In addition, “-omics” experiments 
can estimate the abundances of the many metabolites 
or proteins or regulatory elements at the nodes during 
the steady-state functioning of a cell. However, this in¬ 
formation alone is not sufficient to explain cell function. 

We also need to know the forward and backward rates 
of fluxes LUab and uiba between all nodes a and b, for ex¬ 
ample in metabolic networks Measuring all these 

rates is practically impossible at present, particularly for 
large networks. Third, in structural biology, it is com¬ 
mon to perform computer simulations of the conforma¬ 
tions of biomolecules and infer Markov models among 
metastable states from those simulations ®. Here, com¬ 


puting the populations of the states can be done rapidly, 
whereas computing the kinetic barriers between them is 
much slower. 

In many such problems, the popular approach, espe¬ 
cially for large networks, is to hypothesize a paramet¬ 
ric dynamical model and learn the parameters of this 
model from data. Often, a large amount of data is re¬ 
quired and the parameters learnt are not unique. We 
treat this problem as a matter of inference, in the spirit 
of statistical mechanics, where full distribution functions 
are inferred from a few measured equilibrium-averaged 
properties dmi). We provide a solution employing the 
dynamical analog of the principle of Maximum Entropy, 
Maximum Galiber, seeking the single best model that is 
consistent with under-determined data. 

Mathematically, we seek a Markovian random walker 
on a network that has the maximal path entropy and 
that otherwise satisfies prescribed constraints. Towards 
that goal, we first define a class of walkers that satisfy 
a) a prescribed stationary state distribution {pa} over 
the nodes {a} of a network and b) certain dynamical 
properties defined over the ensemble of stationary state 
paths {r}. Examples of dynamical properties include the 
average distance travelled by the walker per unit time, 
the average number of reactions per per unit time, or the 
average number of amino acid changes per unit time in 
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a constantly evolving protein. 

Below, we first derive the Maximum Caliber Markov 
process. We then illustrate its predictive power with two 
examples, a gene expression network and a network of 
metastable states of a small peptide. 


THEORY 

Consider a Markovian random walker on a directed 
network G with nodes V = {a} and edges E. Assume 
that the random walker has a unique stationary state 
distribution {pa} over nodes {a} that is independent of 
the initial conditions. The instantaneous probability of 
the walker being at a node b at time t, qb{t), is governed 
by 

= X! ^abqa{t) - ^ UJbaqbit) = ^ f2&a<3'a(t)(l) 

a a a 

where the time independent rate of transition ujab from 
node a to node b is non zero if and only (o, b) S E. 

In order to maximize the path entropy, we discretize 
time into time intervals 5t and at a later time, take the 
limit 5t —>■ 0. In the discrete time scenario, the tran¬ 
sition rates uJab with units of inverse time are replaced 
by unitless transition probabilities kab- The matrix k of 
transition probabilities depends on the discretization in¬ 
terval and is given by k = « I -|- ^5t. The second 

approximation is accurate only for <C 1. Here, I is the 
identity matrix. 

Let us define an ensemble {T} of stationary state 
paths r=---—>a—>■&—>- 0 —^d... of the walker 
of total duration T. The entropy of this ensemble is 
St = -J2rPi^)^ogP{r) = -T^pakab^ogkab dEHini)- 
The time normalized path entropy S is 

S = - ^ Pakab'^Ogkab- (2) 

{a,b)GE 

The sum in Eq. is only taken over edges (a, b) £ E oi 
the network. All summations below are also restricted to 
edges (o, b) G E unless otherwise stated. 


The constraints on the paths 

Any discrete time Markov process with {pa} as the 
stationary distribution and kab as transition probabili¬ 
ties satisfies two sets of linear constraints (normalization 
and stationarity). These constraints are understood as 
follows: First, from state i at time t, the system has to 
land at one of the states j at time t + 6t. Second, at 
stationary state, a system in state j at time t + St comes 
from one of the state i. We have, 

^ fcab = 1 V a and '^Pakab = Pb V &. (3) 

b a 


Another important constraint is detailed balance, 
Pakab = Pbkba- Below, we will see how detailed balanced 
constrained can be applied explicitly 

We introduce a node-connectivity variable N such that 
Naa = 0, Nab = 1 */ («, b) G E. The path ensemble 
average of N over all trajectories T is given by dHlISI 

w= E PakabNab- (4) 

ia,b)GE 

(N) is the mean number of transitions in a single time 
step St. Below we use {N) to take the desired limit St ^ 0 
in order to convert the discrete time Markov chain to a 
continuous time Markov process. 

Also, we constrain the ensemble average Kb) of arbi¬ 
trary dynamical rate variables Examples of dynami¬ 
cal constraints include the average distance travelled by a 
particle diffusing on an energy landscape per unit time, 
the number of spin flips per unit time in a spin glass 
model, the average number of reactions per unit time, 
etc. We assume r^a = 0 and = 0 if (o, b) ^ E. = 0 
is a mere convenience and not a strict requirement of our 
development (see supplementary materials for details). 
Given that the network is directed, in general we have 
Nab ^ Nba and ^ 


Maximizing the path entropy, subject to constraints. 


We now maximize the path entropy S (Eq. with re¬ 
spect to transition probabilities kab subject to constraints 
imposed by Eq. and Eq. Using the method of La¬ 
grange multipliers, we write the unconstrained Caliber 

C diini) 


C = S + '^ma i'^Pakab - Pa] + i'^Pakab - 


Pb 


- 7 ^PakabNab “ (^) “ E ^PakabT^b “ Kb) ■ 


a,b 


a,b 


(5) 


We maximize the Caliber with respect to kab and de¬ 
rive a closed form expression for the transition probabil¬ 
ities (see supplementary materials for a detailed deriva¬ 
tion) 

kab = pSt\K^.f^Aab if {a,b) G E and (6) 
V fa9b V Pa 

and kaa = I - J^b^ab- Here pSt = e~^ and A^t, = 
g- Y.i pNab if (^a, b) G E and zero otherwise. Constants fa 
and Pa are determined from self-consistent equations 



( 7 ) 
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Since k = I + ftSt as 6t —>■ 0. We take the limit and 
get 


Aab if (a, b) e E (8) 

and uJaa = -J2b^ab- Eqs. I is the most general result 
of this work. Since constants fa and ga are determined 
from global self-consistent equations Eq. the transition 
rate ujab between any two states a and b depends on the 
structure of the entire network. 


Wa& = M 



Detailed balance 


Application to theories of chemical reaction rates 

Our predicted square-root dynamics has an interest¬ 
ing interpretation in theories of chemical reaction rates. 
In organic chemistry, so-called extra-thermodynamic re¬ 
lationships express empirical observations about how the 
rates and mechanisms of certain types of chemical re¬ 
actions are related to their equilibria. These go by the 
name of the Bronsted relation (often applied to acid-base 
catalysis) (fl^ , the Polanyi relationship for surface catal¬ 
ysis 03), Marcus theory for electron transfer reactions in 
solution (USD, or $-value analysis for protein folding (HZD. 
In short, these approaches express that rates are related 
to equilibrium constants in the form of 


Detailed balance constraint requires Pa^Oab = Pb^ba- 
Thus, if both (a, b) and (6, a) are not in E, ojha = bJab = 0 
and the connection between (a, b) can be removed. Con¬ 
sequently, we assume that the network G is undirected 
i.e. (a, b) G E ^ {b, a) £ E and vice versa. Imposing de¬ 
tailed balance constraints on transition probabilities kab 
is equivalent to constraining symmetrized forms of the 
dynamical variables ^ -|- rl^) (see (jHD for a 

proof). In this case, we have A = A"”" and f = g. Thus 

Wafc = p.f^Aab if (a, b) e E. (9) 

V Pa 

It is easy to check that Eq. satisfies detailed balance. 
Interestingly, the transition rates of a detailed balanced 
Markov process are determined entirely locally from the 
properties of states a and b alone and does not depend 
on the global structure of the network. In contrast, the 
same transition rate depends on the structure of the en¬ 
tire network when detailed balance is not satisfied (see 
Eq.[|). 

When we constrain only the mean jump rate {Aab = 
1 if {a,b) G E and zero otherwise), the transition rates 
are given by the simple expression: 

Wat, = if (a, b) G E. (10) 

V Pa 

This result should be contrasted with other popular func¬ 
tional forms of the transition rates that satisfy a pre¬ 
scribed stationary distribution, e.g. Glauber dynam¬ 
ics dUD and Metropolis dynamics (USD. For example, 
historically, Glauber developed his dynamics to study 
magnetic spins dED. The particular form of the transi¬ 
tion rates was motivated by “a desire for simplicity” (USD. 
The derivation of square root dynamics presented here, 
while not resorting to a ad hoc prescription for transition 
rates, follows the same intuitive principle i.e. finding the 
simplest model that is consistent with a given set station¬ 
ary and dynamical constraints and surprisingly predicts 
a dynamics that is qualitatively different from Glauber 
and Metropolis dynamics. 


uj{x) = cK{x)°‘ 


( 11 ) 


where x is a variable representing a systematic change, 
such as a series of different acids of different pKa's, K{x) 
is the equilibrium constants for that series, and lo{x) are 
the rates of reaction, a is a value that usually ranges 
from 0 to 1, expressing the degree of resemblance of the 
transition state to the products of the reaction. 

Now consider a two state system where state a are re¬ 
actants and state b are products such that pt > Pa- We 
have neglected the population in the transition state since 
it is much smaller than the state populations at equilib¬ 
rium. We have equilibrium constant K = pb/pa > 1- 
If no dynamical constraints are imposed, the transition 
rate uJab of the a ^ h reaction according to the MaxCal 


process is given by uJab = In short, our 

MaxCal approach predicts that a = 1/2 in Eq. 11 is the 
most parsimonious assumption for the reaction mecha¬ 
nism prior to any knowledge of global rate information, 
implying that the transition state is halfway between re¬ 
actants and products. 


ILLUSTRATING EQS. AND [^WITH TWO 
EXAMPLES. 

We test the model predictions on a problem of cor¬ 
related expression of two genes that are transcribed by 
the same promoter and on a problem of the equilibrium 
conformations of a small peptide. 


Expression dynamics of two correlated genes 

It is difficult to infer the underlying regulatory ar¬ 
chitecture of large biochemical networks from stationary 
populations of the components such as mRNAs and pro¬ 
teins. Here, we show how Eq. can accurately predict 
the full chemical master equation (CME) from station¬ 
ary distributions and overall rate parameters. Consider 
a biochemical circuit where two genes are adjacent to a 
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FIG. 1. Schematic of two genes that are simultaneously tran¬ 
scribed by the same promoter. Rates ai, 02 ,and a represent 
the synthesis rates and rates /3i and /32 represent the degra¬ 
dation rates. 


const it utively expressing promoter region. We assume 
that the genes are either transcribed individually or si¬ 
multaneously and that they are degraded individually. In 
this toy example, two stochastic variables are correlated, 
namely the copy numbers of the two mRNA molecules. 
We first construct a CME to mimic the biochemical cir¬ 
cuit in silica (Fig. [^. There are 5 rate parameters in 
the CME: 3 synthesis rates oi, 02 , and a and two degra¬ 
dation rates /3i and j32- If n-i and 712 are the number of 
molecules of the first and the second mRNA, the chemical 
master equation describing the system has the following 
form 

= Q,g ^ 2 ) - p{ni + 1, 712 )) 

at 

+ ai (p(7li, 712 - 1) - p(ni, 772 + 1)) 

-I- a {p{ni - 1,772 - 1) - p(t7i, 772 )) 

-I- /3l ((771 -I- l)p(77l -I- 1 , 772 ) - 77ip(77i,772)) 

+ ^2 ((772 + 1)p(77i,772 + 1) - 772P(77i,772)) 

( 12 ) 

Here, p{ni, 772 ; t) is the instantaneous probability of hav¬ 
ing 77i and 772 molecules of mRNA 1 and 2 respectively 
at time t. The terms correspond to individual synthe¬ 
sis of mRNA 1 and mRNA 2, the simultaneous synthesis 
of both mRNAs, and the degradation of mRNA 1 and 
mRNA 2. 

We assume that we can experimentally estimate the 
joint probability distribution Pss(77i, 772) of the mRNA 
copy numbers at steady state. Additionally, we also as¬ 
sume that we have two reporters that count the total 
number of expression events and the total number of 
degradation events respectively. Note that the reporters 
are agnostic to which of the two RNAs has been syn¬ 
thesized or degraded. With only these three pieces of 
information, can we estimate the transition rate matrix 
of the system? 

We choose the following parameters for the CME: 
( 0 : 1 , 02 , o, ,01, /32) = (Ij 0.5, 2.5, 5,10). The choice ensures 
that the number of any of the two mRNA molecules is 
limited to < 6. This results in a small system size in this 


proof of principle work where the total number of states 
is 6 X 6 = 36. In Fig. [^panels A and B, we show the 
numerically observed joint distribution Pss(77i, 772 ). The 
correlated pattern of expression is apparent in panel B; 
the probability Pss (772 |77i) depends on 771 . The higher the 
value of 77i, the higher 772 values become more probable 
as a result of the correlated expression. 

In the CME, while there are only 5 rate parameters, 
there are 145 possible transitions. There are 30 -I- 30 
transitions that correspond to synthesis of mRNA 1 or 2. 
There are 30 -I- 30 transitions that correspond to degra¬ 
dation of mRNA 1 or 2 and there are 25 transitions that 
correspond to simultaneous synthesis of both mRNAs. 
Each transition has its own rate constant but not all rate 
constants are independent of each other. For example, 
the transition (2, 3) —)■ (2, 2) has a rate constant 3/32 and 
the transition (2,4) —)■ (2,3) has a rate constant 4/^2. 
Moreover, many other transition rates are equal to each 
other, for example, the rates for transitions that increase 
the first mRNA copy number such as (1,2) —(2,2) and 
(3,1) —>■ (4,1) are equal to ai and so on. 

We find that Eq. [^accurately gives the full 145 rate pa¬ 
rameters, without the knowledge of the differential rates 
of synthesis and degradation of the two mRNAs (see sup¬ 
plementary materials for details of the fitting procedure) 
(see Fig. [^. Taken to larger scale, it implies that using 
steady state data on cell-to-cell variability in gene expres¬ 
sion and a few overall kinetic measurements, we can a 
full set of chemical master equations and infer regulatory 
details. This infered model is optimal in the Maximum 
Caliber sense. 


Dynamics of a small peptide 

Second, we study the equilibrium dynamics of 
metastable states of a small peptide comprising 7 alanine 
amino-acid residues. A metastable state is an ensemble of 
geometrically and dynamically proximal microstructures 
that have a significant net population. Classifying pro¬ 
tein structures into their metastable states is an active 
area of research dHHHHSQl). 

As a test, we compare to previous extensive MD simu¬ 
lations of this peptide m, which led to the identification 
of 32 metastable conformations (1^ (see Fig. 1^. From 
that MD simulation, we took the states as defined by 
Jain and Stock and estimated the relative probability Pa 
of each metastable state a as the fraction of time points 
in the full trajectory when the peptide was in that state. 
We also estimated the transition probabilities kab as the 
fraction of the events in the trajectory when the peptide 
was in state a at time t and transitioned into state b at 
time t+6t {5t = 1 ps in a MD simulation of total duration 
T = 800 ns). 

To predict the transition probabilities kab, we first need 
to estimate the transition rate matrix using Eq. [^ In 
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FIG. 2. Panel A: The joint stationary state distribution Pss(ui,n 2 ) when the parameters are set at (ai, a 2 , a,/3i,/32) = 
(1,0.5, 2.5, 5,10). Panel B: The conditional probability p{n 2 \ni) at different values of ni. 
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FIG. 3. Predicted rates vs. toy model values for all 145 
possible transitions in the gene expression model. 
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FIG. 4. Two principal coordinates from molecular dynam¬ 
ics simulations of the alanine-7 peptide m, showing that 
the microscopic structures can be lumped into well separated 
metastable states. The 10 most populated states are labelled 
in decreasing order of state probability. 


order to guess the functional form of the transition rate 
matrix, we need to identify a dynamical constraint. The 
accuracy of our predictions depends on how well the dy¬ 
namical constraint captures the diffusion of the peptide 
on the free energy landscape. Unfortunately, there is no 


systematic procedure to guess a ‘good’ constraint vari¬ 
able. This is a common aspect of Maximum-Entropy 
methods. See (|22H211) for a discussion on the role of con¬ 
straints in maximum entropy methods. 

To guess the dynamical constraint, we make two ob¬ 
servations: the transition rate decreases when 1) when 
the average conformational separation between states is 
increased, keeping the stationary probabilities and the 
free energy barrier constant and 2) when the free energy 
barrier is increased, keeping the average conformational 
distance and the stationary state probabilities constant. 
As a first guess, we only model the effect of geomet¬ 
ric separation and neglect the free energy barrier. We 
choose a simple geometric constraint Vab- for any two 
metastable states a and b and microstates x and y such 

that X G a and y G b, Vab = (/xGa.yefe ^ 

Tab is the mean distance between a randomly chosen mi¬ 
crostate point X in metastable state a and randomly cho¬ 
sen microstate point y in metastable state b. Similarly, 
Xaa is the average distance between any two microstates x 
and y within a macrostate a; it represents the ‘volume’ of 
macrostate a. The distance between any two microstates 
X and y is defined as the Euclidean distance between their 
internal coordinate representation (HOI)- The dynamical 
average {xab) represents the average distance travelled by 
the microstates of the peptide per unit time. 

Using Eq. the numerically estimated stationary 
probabilities, and the chosen geometric constraint, we 
arrive at the functional form of the transition rate ma¬ 
trix r2. The transition matrix had two free parameters, 
/i the time scale and p the Lagrange multiplier associated 
with {xab)- In order to estimate the transition probability 
matrix k from the predicted matrix of transition rates SI, 
we need the discretization time scale a third param¬ 
eter. Since p, and 5t can be combined together, we only 
needed to determine 2 parameters, p and p5t. Note that 
the Maximum Caliber approach guesses the parametric 
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(observed in MD simulation) 

FIG. 5. Panel A: A comparison of transition probabilities {feat} predicted using Eq.j^and those observed in the MD simulation 
for the 7 residue peptide shows good agreement over 5 orders of magnitude. Panel B: Histogram of the logarithm of predicted 
vs observed transition probabilities shows that while the majority of the transition probabilities are predicted very accurately, 
virtually all transition probabilities are predicted within one order of magnitude. 


form of the transition probabilities, one can use any suit¬ 
able numerical technique and experimental information 
to estimate the parameters. Since we have access to mi¬ 
croscopic data, in this proof of principle study, we used 
the entire transition probability matrix to fit the two free 
parameters. Using multiple simulated annealing runs to 
minimize the total error between known and predicted 
transition probabilities, we found that the best fits were 
at fiSt = 45 ± 5 and p = 4.8 ± 0.6. Fig. shows that 
the rates predicted by the Max Cal approach quite accu¬ 
rately capture the correct values obtained from the full 
MD simulation, over 5 orders of magnitudes of rates. 


DISCUSSION 

Here, we describe theory that takes a given network, 
steady-state populations on its nodes, and a couple of 
global dynamical constraints, and finds the microscopic 
transition rates among all the nodes. We do this us¬ 
ing Maximum Caliber, a Maximum-Entropy-like princi¬ 
ple for dynamical systems. A main finding is that the 
MaxCal transition rates are proportional to the square 
root of ratios of the state populations. We illustrate 
our results on a toy gene expression network and pep¬ 
tide conformations, for which we know the correct rates 
in advance. We believe this treatment could be useful in 
many areas of network modeling, including in spin-glass 
models of the immune system (H USD, colloidal assem¬ 
blies ®, neuronal networks m, master-equation models 
of noisy gene expression (j24l |26]), and in the browsing 
behavior of web crawlers on the internet (EZl). 
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SUPPLEMENTARY MATERIALS 
Deriving the maximum Caliber transition rates 


where ^ = 6"*““^, A;, = and Wat, = 6“'*'^“'’“^’'“'’. 

For a given value of 7 and p, the Lagrange multipli¬ 
ers /3s and As are determined by self consistently solving 
Eqs. We have 

^fca6=l^^Waf,At, = §^ 
b b 

and ^ Pakab =Pb^'^ '^ba/3b = ^ • (14) 


To further simplify Eqs. 14 let usidentify A and /3 as the 
vectors of Lagrange multipliers and define 'D[x]a = a 
non-linear operator on vectors x. We have 


WA = V[P] and W'^/3 = V[X] 

=7 I?[WA] = ^ and I?[W^^] = A (15) 


Now, we write W = I + p6tA where e“'>' = p6t, Aab = 
0 if (a, 6) ^ E, and Aab = when (a, 5) S E. We 

have 

= Pa and = Aa (16) 

Aa + pot fa Pa + POtpa 

where 


/ = AA and g = A'^/3. (17) 


We solve these two equations algebraicly recognizing that 
fa and Pa do not directly depend on Aa and j3a- The 
positive roots are 


Aa 


faPaPp + 4pa - StfaPaP 
‘^9a 

VJay/^VSt^faGaP^ + ^Pa - StfaPaP 

2fa 


(18) 

(19) 


Given that we’re only interested in transition probabili¬ 
ties kab up to orders of 6t and kab = pSt^XbAab already 
has a St term, we only need to the zeroth order terms for 


simultaneously for Aa and /3a in terms of fa, Pa, Pa, and 
dt. We then take only the zeroth order terms in St (see 
supplementary materials for details). We have 


A and (3 as St ^ 0. We solve the algebraic equations 16 


A 


a 



and /3a 



( 20 ) 


Eq. 20 and Eq. 17 can be self-consistently solved for A 
and f3 and we can obtain / = AA and g = A'^p. 


For notational simplicity, we derive Eq. for one dy¬ 
namical constraints Vab- Differentiating the Caliber in 
Eq.0 with respect to kab setting the derivative to 
zero, 

PaO^Ogkab + 1) = maPa + hPa “ pNabPa “ PTabPa 

^ kab = —XbWab (13) 

Pa 


When dynamical constraints are such that Vaa f 0 

As mentioned in the main text, we have assumed for 
convenience that the dynamical constraints Xab is such 
that Taa = 0. Here, we show that when Xaa 7 ^ 0, the 
maximum entropy problem is equivalent to constraining 
a modified constraints = Tab — \ ijaa + ?'&b)- 

































We start with recognizing W = J + fiStA. as above 
where J is a diagonal matrix such that Jaa = . We 

have 


Pa 


jaK + pStfa 


= Pa and 


Pa 


jaPa + pStga 


= Aa (21) 


where ja = Jaa and / = AA and g = A'^/3. Again 
solving to zeroth order in St, we get 


A„ = 




QclJcl 


— and Pa = 


Pa9a 

faja 


( 22 ) 


2. Synthesis or degradation of mRNA 2: ni = n\ and 

712 = n.2 ± 1 

3. Simultaneous synthesis: n\ = ni+1 and rij = n 2 +l 

As mentioned in the main text, the ‘experiments’ tell 
us the total number of degradation and synthesis events 
per unit time but do not have the ability to distinguish 
between the two mRNAs. We constrain two quantities, 
the number of degradation events per unit time and the 
number of synthesis events per unit time. Accordignly, 
we set A in Eq. as follows 


Finally, the transition probability kab is given by 


kab = pSt^XbC-P^^'- = 

Pa V V JaPb y JaJb 

\ Pa \ fagb 


(23) 


where i 


' ab 


= Tab - I (r 


■Tbb)- Comparing Eq. 


23 


to 


Eq. it is clear that the problem of constraining a dy¬ 
namical constraint Tab such that Vaa is non-zero is equiv¬ 
alent to constraining a modified dynamical constraint 
^Ib = Pab-^ (raa + nb)- Note that by definition, rl^ = 0. 
Thus, for convenience, we assume that this transforma¬ 
tion is already performed and raa = 0. 


Constructing and fitting the maximum entropy 
Markov proces for the gene network 

From the numerical experiments, we obtain accurate 
estimate of the stationary state probability distribution 
Pss{ni,n 2 ). Let a = (ni, 712 ) and b = ( 77 }, 772 ) be any two 
states of the system. We know that there is a directed 
edge from state a to state b iff 

1. Synthesis or degradation of mRNA 1: 771 = tt-I ± 1 
and 772 = n 2 


1. No edge between nodes a and b: Aab = 0 if (a, b) ^ 
E 

2. Synthesis events: Aab = 77 if 77 i = 77 J and 772 -I- 1 = 

772 or 77i -I- 1 = 772 = 772 or Til + 1 = 

77 I and 772 -f 1 = 772 

3. Degradation of mRNA 1: Aab = tt-iC if ni = 77^ -I- 
1 and 772 = 77-2 

4. Degradation of mRNA2: Aab = 772 C if 77i = 

77^ and 772 = 772 -|- 1 


Here, 77 > 0 and C ^ 0 are exponentials of Lagrange 
multipliers similar to those used in the main text. The 
Lagrange rate constant g, in Eq. is assumed to be ab¬ 
sorbed in the Lagrange multipliers. The next step is to 
determine /, g, A, and p using numerically estimated 
Pss( 711 , 772 ) and A. In order to determine the transition 


rate matrix El for any value of 77 and we solve Eq. 17 
and Eq. self consistently. 

Given that we have access to all rate constants in this 
proof of principles work, we minimize the error between 
the known rate constants of Eq. 12 and those predicted 
by Eq. iby varying 77 and (f using a simulated annealing 
protocol. We find that 77 « 0.4 and C ~ 26.5 resulted 
in the best agreement between the known and the pre¬ 
dicted rate constants. Multiple simulated annealing runs 
predicted rates that were identical to the ones shown in 
the main text. Note that although we used the entire 
transition matrix to learn the Lagrange multipliers, it 
is equally possible to learn them from global dynamical 
constraints. 

















